"""
https://vcfw.org/pdf/Department/Physics/Fourier_series_python_code.pdf
# scipy.signal.square (x, duty=0.5)
# scipy.signal.sawtooth(x, width=1)
# scipy.signal.triang(M, sym=True)
"""
import numpy as np
from scipy.signal import square, sawtooth, triang
import matplotlib.pyplot as plt
from scipy.integrate import simps


def a0(L, x, y):
    return 2. / L * simps(y, x)


def an(L, x, y, n):
    return 2.0 / L * simps(y * np.cos(2. * np.pi * n * x / L), x)


def bn(L, x, y, n):
    return 2.0 / L * simps(y * np.sin(2. * np.pi * n * x / L), x)


# Fourier series analysis for a square wave function
L = 4  # Periodicity of the periodic function f(x)
freq = 4  # No of waves in time period L
duty_cycle = 0.5
samples = 1000
terms = 100
# Generation of square wave
x = np.linspace(0, L, samples, endpoint=False)
y = square(2.0 * np.pi * x * freq / L, duty=duty_cycle)
# Calculation of Fourier coefficients
a0 = 2. / L * simps(y, x)
an = lambda n: 2.0 / L * simps(y * np.cos(2. * np.pi * n * x / L), x)
bn = lambda n: 2.0 / L * simps(y * np.sin(2. * np.pi * n * x / L), x)
# sum of the series
s = a0 / 2. + sum([an(k) * np.cos(2. * np.pi * k * x / L) + bn(k) * np.sin(2. * np.pi *
                                                                           k * x / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(x, s, label="Fourier series")
plt.plot(x, y, label="Original square wave")
plt.xlabel("$x$")
plt.ylabel("$y=f(x)$")
plt.legend(loc='best', prop={'size': 10})
plt.title("Square wave signal analysis by Fourier series")
plt.savefig("fs_square.png")
plt.show()

# Fourier series analysis for a sawtooth wave function
L = 1  # Periodicity of the periodic function f(x)
freq = 2  # No of waves in time period L
width_range = 1
samples = 1000
terms = 50
# Generation of Sawtooth function

x = np.linspace(0, L, samples, endpoint=False)
y = sawtooth(2.0 * np.pi * x * freq / L, width=width_range)
# Calculation of Co-efficients

an = lambda n: 2.0 / L * simps(y * np.cos(2. * np.pi * n * x / L), x)
bn = lambda n: 2.0 / L * simps(y * np.sin(2. * np.pi * n * x / L), x)
# Sum of the series
s = a0 / 2. + sum([an(k) * np.cos(2. * np.pi * k * x / L) + bn(k) * np.sin(2. * np.pi *
                                                                           k * x / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(x, s, label="Fourier series")
plt.plot(x, y, label="Original sawtooth wave")
plt.xlabel("$x$")
plt.ylabel("$y=f(x)$")
plt.legend(loc='best', prop={'size': 10})
plt.title("Sawtooth wave signal analysis by Fouries series")
plt.savefig("fs_sawtooth.png")
plt.show()

# Fourier series analysis for a Triangular wave function
L = 1  # Periodicity of the periodic function f(x)
samples = 501
terms = 50
# Generation of Triangular wave
x = np.linspace(0, L, samples, endpoint=False)
y = triang(samples)
# Fourier Coefficients
a0 = 2. / L * simps(y, x)
an = lambda n: 2.0 / L * simps(y * np.cos(2. * np.pi * n * x / L), x)
bn = lambda n: 2.0 / L * simps(y * np.sin(2. * np.pi * n * x / L), x)
# Series sum
s = a0 / 2. + sum([an(k) * np.cos(2. * np.pi * k * x / L) + bn(k) * np.sin(2. * np.pi *
                                                                           k * x / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(x, s, label="Fourier series")
plt.plot(x, y, label="Original Triangular wave")
plt.xlabel("$x$")
plt.ylabel("$y=f(x)$")
plt.legend(loc='best', prop={'size': 10})
plt.title("Triangular wave signal analysis by Fouries series")
plt.savefig("fs_triangular.png")
plt.show()

# Fourier series analysis for a sawtooth wave function
# User defined function
L = 1.0  # half wavelength, Wavelength=2L
freq = 2  # frequency
samples = 1001
terms = 300
# Defining sawtooth function
x = np.linspace(-L, L, samples, endpoint=False)
f = lambda x: (freq * x % (2 * L) - L) / L
# Fouriers coefficients

a0 = 1. / L * simps(f(x), x)
an = lambda n: 1.0 / L * simps(f(x) * np.cos(1. * np.pi * n * x / L), x)
bn = lambda n: 1.0 / L * simps(f(x) * np.sin(1. * np.pi * n * x / L), x)
# Series sum
xp = 4 * x
s = a0 / 2. + sum([an(k) * np.cos(1. * np.pi * k * xp / L) + bn(k) * np.sin(1. * np.pi
                                                                            * k * xp / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(xp, s, label="Fourier series")
plt.plot(xp, f(xp), label="Original sawtooth wave")
plt.legend(loc='best', prop={'size': 10})
plt.savefig("saw_ud.png")
plt.show()

# Fourier series analysis for a square wave function
# User defined function
L = 1.0  # half wavelength, Wavelength=2L
freq = 2  # frequency
samples = 1001
terms = 300
# Generating Square wave
x = np.linspace(-L, L, samples, endpoint=False)
F = lambda x: np.array([-1 if -L <= u < 0 else 1 for u in x])
f = lambda x: F(freq * x % (2 * L) - L)
# Fourier Coefficients
a0 = 1. / L * simps(f(x), x)
an = lambda n: 1.0 / L * simps(f(x) * np.cos(1. * np.pi * n * x / L), x)
bn = lambda n: 1.0 / L * simps(f(x) * np.sin(1. * np.pi * n * x / L), x)
# Series sum
xp = 4 * x
s = a0 / 2. + sum([an(k) * np.cos(1. * np.pi * k * xp / L) + bn(k) * np.sin(1. * np.pi
                                                                            * k * xp / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(xp, s, label="Fourier series")
plt.plot(xp, f(xp), label="Original Square wave")
plt.legend(loc='best', prop={'size': 10})
plt.savefig("square_ud.png")
plt.show()

# Fourier series analysis for a Arbitrary waves function
# User defined function
L = 1.0  # half wavelength, Wavelength=2L
freq = 2  # frequency
samples = 1001
terms = 300
# Generating wave
x = np.linspace(-L, L, samples, endpoint=False)
F = lambda x: np.array([u ** 2 if -L <= u < 0 else 1 if 0 < u < 0.5 else 0
                        for u in x])
# F=lambda x: abs(np.sin(2*np.pi*x))
f = lambda x: F(freq * x % (2 * L) - L)
# Fourier Coefficients
a0 = 1. / L * simps(f(x), x)
an = lambda n: 1.0 / L * simps(f(x) * np.cos(1. * np.pi * n * x / L), x)
bn = lambda n: 1.0 / L * simps(f(x) * np.sin(1. * np.pi * n * x / L), x)
# Series sum
xp = x
s = a0 / 2. + sum([an(k) * np.cos(1. * np.pi * k * xp / L) + bn(k) * np.sin(1. * np.pi
                                                                            * k * xp / L) for k in range(1, terms + 1)])
# Plotting
plt.plot(xp, s, label="Fourier series")
plt.plot(xp, f(xp), label="Original wave")
plt.legend(loc='best', prop={'size': 10})
plt.savefig("arb_ud.png")
plt.show()
